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Charge-Balanced Minimum-Power Controls for 
Spiking Neuron Oscillators 

Isuru Dasanayake and Jr-Shin Li 



Abstract 

In this paper, we study the optimal control of phase models for spiking neuron oscillators. We 
focus on the design of minimum-power current stimuli that elicit spikes in neurons at desired times. We 
furthermore take the charge-balanced constraint into account because in practice undesirable side effects 
may occur due to the accumulation of electric charge resulting from external stimuli. Charge-balanced 
minimum-power controls are derived for a general phase model using the maximum principle, where the 
cases with unbounded and bounded control amplitude are examined. The latter is of practical importance 
since phase models are more accurate for weak forcing. The developed optimal control strategies are 
then applied to both mathematically ideal and experimentally observed phase models to demonstrate 
their applicability, including the phase model for the widely studied Hodgkin-Huxley equations. 

Index Terms 

Spiking neurons, Phase models, Optimal control, maximum principle, pseudospectral method. 

I. Introduction 

The electrical activity of a nervous system and its ability to respond to external electrical 
signals have been long-standing subjects of active research. The resulting insights have led to 
the innovation of therapeutic procedures for a wide variety of neurological disorders. Deep brain 
stimulation is one such method applying electrical pulses to inhibit pathological synchrony among 
the neurons 0] and is clinically approved in many countries for the treatment of Parkinson's 
disease, essential tremor, and Dystonia 0, 0. The cardiac pace maker is another example in 

Isuru Dasanayake and Ir-Shin Li are with the Department of Electrical and Systems Engineering, Washington University, St. 
Louis, MO, 63130, USA. e-mail: dasanayakei@seas.wustl.edu, jsli@seas.wustl.edu 
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medical practices that employs electric pulses to stimulate nervous tissues in order to regulate 
a patient's heart rate JH, JH. In these and many other neurological applications, the use of 
low power electrical stimuli is desired because, for example, high power stimuli are harmful to 
biological tissues and the reduction of power consumption in a neurological implant is essential 
in order to reduce its sizes and lengthen its lifetime. In addition, it is of clinical importance to 
ensure that any external inputs, e.g., currents, applied to stimulate neurons are charge-balanced. 
That is, the net amount of the electric charge injected into a neuron over one oscillation cycle 
should be kept zero, because high levels of the charge accumulation may trigger irreversible 
electro-chemical reactions, resulting in damage of neural tissues and corrosion of electrodes J6]|. 

Many mathematical models have been developed to capture the periodic activities of neuron 
oscillators 0, flU, 0, ifTOl and a well established example is the phase response curve (PRC), 
which quantifies the asymptotic phase shift of an oscillator due to an infinitesimal perturbation of 
its state [1 1 J. A phase model accurately approximates the behavior of the corresponding full state- 
space system in the neighborhood of its periodic orbit lfT2l . Due to their simplicity, phase models 
are very popular for modeling and analyzing the dynamics of neuron oscillators. For example, 
the patterns of synchrony resulting from the dynamics of an arbitrary network of oscillators 
with weak coupling were analyzed using phase models [TT3l . Ifl4l . and a chain of coupled phase 
oscillators has been used to model the lamprey spinal generator for locomotion [fT5l . In these 
studies, the inputs to the oscillatory systems were initially defined, and the dynamical responses 
of neuron populations were analyzed in detail. Recently, as an alternative objective, control 
and dynamical systems approaches have been used to manipulate neural activities in a desired 
way. For instance, minimum-power controls for spiking neurons at specified time instances were 
derived for some mathematically ideal phase models lfT6ll . ifTTl and charge-balanced controls 
were calculated using a numerical shooting method lfT8ll . Controllability of a network of neurons 
described by phase models has also been investigated lfT9l . 

In this article, we consider a general phase model and derive charge-balanced minimum-power 
controls for spiking a neuron oscillator at a desired time instance different from its natural spiking 
time. Both cases of unbounded and bounded control amplitude are examined. The latter is of 
fundamental and practical importance since there exist physical limitations on medical equipment 
and safety margins for neural tissues and, more importantly, because phase models are valid under 
weak forcing. We show that the bounded optimal control has switching characteristics synthesized 
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by the unbounded optimal control and the given control bound. The developed optimal control 
strategies are then applied to both mathematically ideal, such as sinusoidal, and experimentally 
observed, such as Hodgkin-Huxley, PRC's to demonstrate their applicability. In addition, we 
characterize the range of possible spiking times with respect to the given control amplitude for 
several phase models. Moreover, we apply the optimal controls derived from the reduced phase 
model to the corresponding full state-space model to verify the consistency of these models 
through the reduction and the robustness of our optimal control techniques. Such an important 
validation is missing in the literature. 

This paper is organized as follows. In Section [HI we present the optimal control problem of 
spiking a general phase oscillator. We find the charge-balanced minimum-power control for a 
prescribed spiking time with and without a constraint on the control amplitude by using the 
maximum principle. In Section Unl we apply the derived optimal control strategies to several 
commonly-used phase models and present the optimal solutions and numerical simulations. In 
particular, we calculate optimal controls for experimentally observed PRC's including Morris- 
Lecar and Hodgkin-Huxley PRCs. These optimal controls produced by the maximum principle 
are verified by the Legendre pseudospectral computational method ll20l . 

II. Optimal charge-balanced controls for spiking neurons 

In systems theory, a nonlinear oscillator is described by a set of ordinary differential equations 
that has a stable periodic orbit. This system of equations can be reduced to a single first order 
differential equation, which is valid while the state of the full system remains in a neighborhood 
of its unforced periodic orbit |fTT|. This reduction allows us to represent the dynamics of a 
weakly forced oscillator by a single phase variable that defines the evolution of the oscillation. 
Consider a time-invariant system x = f(x, I), where x(t) G M n is the state and I(t) G R is the 
control, which has an unforced stable attractive periodic orbit 7(2) = j(t + T) homeomorphic 
to a circle, satisfying 7 = f(j, 0). We can represent this system in a phase-reduced form as 

9 = f(9) + g(9)I(t), (1) 

where 9 is the phase variable, / and g are real- valued functions, and I(£) 6 I is the control 
ifTTfl. |fT2|. One complete oscillation of the system corresponds to 9 G [0, 2tt). The function / 
gives system's baseline dynamics and g is known as the phase response curve (PRC), which 
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describes the infinitesimal sensitivity of the phase to an external control input. In the case of 
neural oscillators, I represents an external current stimulus and / is referred to the instantaneous 
oscillation frequency in the absence of any external input, i.e., 1 = 0. Neuron spiking occurs 
when the oscillator evolves through one complete cycle. As a convention, the occurrence of spikes 
takes place at 9 = 2nn, where n = 0,1,2, .. .. We consider spiking a neuron at a prescribed time 
with a minimum-power stimulus and, furthermore, intend to find a charge-balanced one in order 
to minimize the side-effects cause by the accumulation of electric charge. The design of such 
charge-balanced minimum-power current stimuli for spiking neurons gives rise to a constrained 
optimal steering problem for a single-input nonlinear system of the form 

min I I(t) 2 dt, 

m Jo 

s.t. 6 = f(6)+g(6)I(t), 
P = I(t), 

9(0) = 0, 6{T) = 2ix, (P) 

p(0) = 0, p(T) = 0, 

\I(t)\<M, 

where M G M + defines the bound of the control amplitude, and the time-dependent variable 
p(t) = f* I (a) da, with boundary conditions p(0) = p(T) = 0, is introduced to accommodate the 
charge-balanced constraint. In the following, we first consider the case of an unbounded control, 
namely, M = oo, and then extend the result to the case when the control is bounded. 

A. Charge-Balanced Minimum-Power Control with Unconstrained Amplitude 

Relaxing the amplitude constraint by letting M = oo, we apply the maximum principle to 
characterize the extremal trajectories. The Hamiltonian of the optimal control problem (0 is 
given by 

H = X I 2 + X(f(e)+g(e)I)+fj,I, (2) 

where A , A, and p are Lagrange multipliers associated with the Lagrangian, system dynamics, 
and the charge-balanced constraint, respectively. Here we consider normal extremals which are 
found by taking A ^ 0. Note that more specific abnormal extremals found by letting A = can 
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be analyzed according to the expressions and properties of the functions / and g. Our derivations 
here are made for the general phase model, and therefore these extraordinary cases are omitted. 
Abnormal extremals are in general uncommon in phase models, and none of the phase models 
considered in this paper has an abnormal extremal (see Remark [2] in Section UlI-AI) . Therefore, 
without loss of generality, we let A = 1. The optimality condition from the maximum principle 
demands that §7 = along the optimal trajectory, which yields 

/ = _*»W±£. (3 ) 

The adjoint variables A and \i are solutions to the time- varying differential equations A = — ^ 
and fi — — Together with © these equations can be written as 

• df{0) X(Xg + fi) dg(6) 

A =0, (5) 

which implies that fj, is a constant. In addition, since the Hamiltonian is not explicitly dependent 
on time, H is a constant along the optimal trajectory. Hence, we let H = c, V< G [0,T]. This 
can also be seen from the transversality condition of the maximum principle. 

It follows that the optimal multiplier A can be found from © by substituting © for /. Then, 
solving for A yields 

> _ -M + 2f±2y/P-g»f-g*c 

g 2 

Here we will choose the negative square root because the positive case corresponds to a backward 
evaluation of the phase, which would invalidate the phase model. The phase velocity equation 
along the optimal trajectory can then be found by using ©, ©, and (OQ), resulting in 

o = VP - g»f - g 2 c (7) 

In addition, substituting © into © gives rise to the optimal control I* in terms of the two 
constants [i and c, 

r = -/ + V f 2 - gyf - g 2 c 
g 

For a given spiking time T, the constants c and /i can be determined from © by separation of 
variables together with the charge-balanced constraint written as J 2w ^-j^-dO = 0, which yields 

r 2n i 

T= / =d9, (9) 

Jo \/f 2 -g^f-g 2 c 
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and 

*-f+/r-9rf-9>c„ Q (10) 
9Vf 2 -9f*f ~9 2 c 

Now the optimal control is completely classified by ©, because the constants /j and c can be 
found from ([9]) and (flOl) for any specified spiking time T. 

Remark 1: In the absence of the charge-balanced constraint, corresponding to // = 0, it is 
sufficient to characterize the optimal control by © and ©. 

5. Charge-Balanced Minimum-Power Control with Constrained Amplitude 

In practice, the feasible amplitude of the stimulus is limited, and phase models are valid only 
for weak forcing. Therefore, spiking neurons with controls of bounded amplitude is of practical 
importance. In this case where we assume that |/| < M, V t G [0, T], the minimum and maximum 
possible spiking times for a neuron system can be determined. It is easy to see that for a given 
bound M > 0, the minimum spiking time is achieved by 

f M, g(6) > 
It ={ (ID 

■* mm 1 

{ -M, g{9) < 0, 

which keeps the phase velocity at its maximum. The minimum spiking time for a given value 
of M, denoted by T^ in , is then given by 

T ™" = / HO) + g mM M + I M -MM*' (12) 

eeA e&B 
where the sets A and B are defined as 

A = {6\ g(0) > 0,0 < 6 < 2tt}, 
B ={9\ g{9) < 0,0 < 9 < 2tt} . 

Symmetric to the minimum spiking time, the maximum spiking time, denoted T^ ax , for the bound 
M can be found by applying the opposite control —I^ m _ n , for M < min{|^y| : 9 G [0, 2tt)}, 
and it is given by X^ ax = T m ^. Note that arbitrarily large spiking times are achievable if the 
bound M > min{|gg| : 9 G [0,2tt)}. 

It is obvious that if < M, V 9 G [0, 2ir), then the amplitude constraint is inactive and I* 

as in © is the charge-balanced minimum-power control. While |/*| > M for some 9 G [0, 2%), 
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it is sufficient to consider the case when I* > M because the case I* < —M is symmetric. 
Suppose that I* > M for 9 G (9x, 9 2 ) C [0, 2n), we now show that the bang control I — M is 
optimal for # G [#i , 2 ] . Since the Hamiltonian © is a convex function of /, / = M is then the 
minimizer when I* > M for 9 G #2]- In this case, we have, from ©, the Lagrange multiplier 

f(9) + M<?(0) ' 1 " } 

which satisfies the adjoint equation (HJ), and hence 1(9) = M is optimal for 9 G [^i, 62]. Similarly, 
the same approach can be used to show that I = — M is optimal on the interval over which I* < 
—M. Therefore, the charge-balanced minimum-power control with limited control amplitude M 
is of the form with switching characteristic 

' — M, I*(9) < -M 

I*m{6) = \ I*(9), -M <I*{9) <M (14) 

M, 1(9)* > M. 

The switching phases 9 G [0, 2tt) such that I* (9) = -M or I* (9) = M can be computed (see the 
examples in Section [ill]) and the required parameter values \i and c can be calculated according 
to the specified spiking time and the charge-balanced constraint from the equations 

T =Cw)TmF M M <15) 



and 

f2n j 



M 



f(9)+g(9)P M 
III. Example 



d9. (16) 



We now apply our optimal control strategies to several commonly-used phase models char- 
acterized by various PRC's, including mathematically ideal models, such as sinusoidal PRC, 
SNIPER PRC, and theta neuron PRC, as well as more realistic phase models such as Hodgkin- 
Huxley and Morris-Lecar PRC's. These mathematically ideal phase models are approximations to 
full state-space models at certain bifurcation points, whereas Hodgkin-Huxley and Morris-Lecar 
phase models are obtained numerically by perturbing their periodic orbits using unit impulses. 
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A. Sinusoidal Phase Model 



The sinusoidal phase model is characterized by a sinusoidal PRC ifTTTl . 



9 = cu + Zd(sm8)I, 



(17) 



where uo is the natural oscillation frequency of the system, z d is a model-dependent constant, 
and / is the external stimulus. This is a type II PRC, with both positive and negative regions, 
which results from a periodic orbit near the supercritical Hopf bifurcation IfTTTl . and occurs in 
neuron models such as the abstracted FitzHugh-Nagumo neuron model [21 J. Neurons described 
by this phase model spike periodically with the natural period T = 2tt/lu in the absence of any 
external input. 

Observe from © that with / and g as defined above, I* (9) is anti- symmetric around 9 = ir, 
namely, I* (9) = —I* (9 + it) for < 9 < it. Therefore, the charge-balanced constraint is 
automatically fulfilled for the sinusoidal phase model. As a result, we let fx = 0. 

1 ) Unbounded Control for Sinusoidal Phase Model: Substituting f = u> and g = Zd sin 9 with 
[i = into ® and ©, the optimal control for spiking a sinusoidal neuron at time T is 



A simple example is used to demonstrate these results. For a neuron with the natural oscillation 
frequency cu = 1 and Zd — 1, the optimal controls for the desired spiking times T = 4 and 
T = 9, smaller and greater, respectively, than the natural spiking time T = 2n are shown in 



Remark 2: Abnormal extremals in general do not exist in phase models. Consider the case 
of abnormal extremals for the sinusoidal phase model, where the multiplier A = 0. Then, the 
Hamiltonian as in (O is given by H = \cu + \zd(sin9)I + /if, and the optimality condition of 
the maximum principle gives 




(18) 




(19) 



Fig. |l(a)| The corresponding optimal phase trajectories are depicted in Fig. |l(b)[ 



-— = \z d sin 9 + /i = 0. 



(20) 



Differentiating this equation with respect to time, we obtain 



Xzdicos 9)9 + \z d sin 9 + fi = 0. 



(21) 
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Time (t) 

(a) 



■■■'1=0 
1=4 




Time (t) 

(b) 



i ■) 
- 1 -i 



Fig. 1. Sinusoidal phase model with u — 1 and za — 1. |(a)| Unbounded charge-balanced minimum-power controls for the 
spiking times T = 4 and T = 9, |(b)| Optimal phase trajectories following the optimal controls. 



Substituting (fTTT) . ©, and © into (|2"TI) for 0, A, and /i, respectively, yields 

cuXzdCosO = 0. (22) 



Abnormal extremals must satisfy (|22|) . and it is clear that (122)) holds only when A = 0. This 
leads to ji = from (|20l) . which, together with A = 0, violate the nontriviality condition of the 
maximum principle. 

2) Bounded Control for Sinusoidal Phase Model: As presented in Section III-BL with the 
amplitude constraint \I(t)\ < M, Vt E [0,T], there exists a range of times at which a neuron 
can be fired. According to (fT2l) for z d > 0, the minimum possible spiking time is given by 

I 4 tan" 1 {z d M/^-z 2 d M 2 +w 2 } 



A min — ^ 'M 



Observe from (TTTI) that when M > u/zd, arbitrarily large spiking times can be achieved by 
making 6 arbitrary close to zero. Therefore, the maximum spiking time T^ ax is given by T~*f 
for M < co/zd and the value of T^ ax is infinity for M > ui/z d . It follows that the assignment of 
the spiking time to any T E K„i^!L] 1S feasible with the control I* M as in CHI). Obviously, 
if the amplitude of the unbounded optimal control satisfies |/*| < M for all t E [0, T], or 
equivalently, V# E [0, 2n), then the amplitude constraint is inactive and I* will be the charge- 
balanced minimum-power control for this bound M. There exists a shortest possible spiking 
time achievable by /* under the bound M, namely (see Appendix |A]) 

T^n= f 1 = -Ad. (23) 

'u 2 + z d M(z d M + 2u) sin 2 6 
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The maximum spiking time achieved by I* is T^ ax = T^ in \ M= _ M for M < ui/z d and T^ ax = oo 
for M > u/z d (see Appendix®. Note that T^ n < T^ in < T% ax < T™ ax and a spiking time 
T G (0,X^ n ) U (T^ ax , oo) cannot be achieved with the bound M. In order to properly classify 
the feasible spiking ranges and associated controls, we consider the two cases, where M < co/zd 
and M > co/zd- 

Case I: (M < uj/z d ) For a desired spiking time T G [T*{ n ,T^*J, the charge-balanced 
minimum-power control I* M is characterized, according to (fl~4~l) . by switching between I* and 

M, 

' i* o < e < d x 

M 9 l <9<9 2 

r e 2 < e < e 3 (24) 

—M 9 3 <9<9 4 
I* 8 4 <6< 2tt, 

in which 9 X = sm~ 1 [-2Mu/(z d M 2 + z d c)\, 9 2 = 7r - X , 3 = 7T + 6>i, and 6» 4 = 2vr - 6»i (see 
Appendix |A)) . The constant c can be computed according to the desired spiking time T, as in 
(fl"5T) . through the relation 



T = / =d9 + / -d^. (25) 

Jo ^u 2 - cz 2 sin 2 9 J 0l uj + z d M sin 9 

The spiking time T G [T^,^] can be optimally achieved by the control J*, and for T 6 
[T^T^j] the optimal control is given by substituting M = —M in the expressions (1241) and 
(1251) . i.e., I* M . A summary of the optimal (minimum-power) spiking scenarios for a prescribed 
spiking time of a neuron governed by the sinusoidal phase model (flTI) is illustrated in Fig. [2] 
Fig. [3] shows the optimal controls for spiking a sinusoidal neuron with to = 1 and ^ = 1 at 
T = 4.7, 5, 8, 10 with the control bound M = 0.6 < u;/^ = 1. These spiking times are chosen 
to cover all possible spiking scenarios depicted in Fig. [2] For this particular example, we select 
the cases of both T <T = 2n/u and T > T , where T = 4.7 G [T^ n , T^J, 5, 8 G [T^ n , T£ x ] 
and 10 G [T^T^J. 

Case //: (M > co/zd) In this case, arbitrarily large spiking times are possible because the 
system can be driven arbitrarily close to the equilibrium point 9 = 0. Analogous to the previous 
case, if the desired spiking time is T G [J^ n , T£ in \, then the switching control I* M , as given in 
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Not feasible Not feasible 

with bound M I M I I_ M with bound M 



T M J 1 ' rpf J,M 

min mm -'max max 

Spiking Time (T) 

Fig. 2. A summary of the optimal control strategies for the sinusoidal PRC model for M < zj/oj 

0.8 

0.6 
0.4 
| 0.2 

I 
E 

§■-0.2 
-0.4 
-0.6 
-0.8 ( 

Fig. 3. Optimal bounded controls with bound M — 0.6 for sinusoidal phase model (with u = 1, Za = 1) to elicit spikes at 
T = 4.7, 5, 8, 10. 

([24]) . will be optimal, and for T G [T^* n , oo) the control I* will be optimal. A summary of optimal 
(minimum-power) spiking scenarios for this case is illustrated in Fig. HI and Fig. [5] shows the 
optimal controls for spiking a sinusoidal neuron with tu = 1 and Zd = 1 at T = 3.5, 4, 8, 12 given 
the control bound M = 1.5 > to/z d = 1. As in the previous case, these spiking times are chosen 
to cover all possible spiking scenarios depicted in Fig. 0] for example, T = 3.5 G [T^f n , T m ;J, 
and 4, 8, 12 G [T^- n , oo]. Note that in this case T^ ax = T^ ax = oo. 

B. SNIPER Phase Model 

SNIPER phase model is characterized by f(6) = to and the PRC g(6) = z d (l - cos 9) IfTTH . 
This phase model is derived from a SNIPER bifurcation (saddle-node bifurcation of a fixed 
point on a periodic orbit) which can be found on type I neurons ll22l like the Hindmarsh-Rose 
model. The charge-balanced minimum-power control for unbounded control amplitude can be 
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Fig. 4. A summary of the optimal control strategies for the sinusoidal PRC model for M > u/zd 
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Fig. 5. Optimal bounded controls with bound M = 1.5 for sinusoidal phase model (with u = 1, Zd = 1) to elicit spikes at 
T = 3.5, 4, 8, 12. 



readily calculated according to ®, ©, and (flOl) . and for bounded control amplitude the control 
is calculated according to (fl4l) . (TT5T) . and (TT6l) using these / and g functions. In Fig. |6(a)| and 
\6(b)\ we show unbounded optimal controls in the absence and presence of the charge-balanced 
constraint and the resulting trajectories of a SNIPER neuron with u — 1 and za — 1. Note that 
the optimal controls without considering the charge-balanced constraint are obtained by taking 
fx = 0. Fig. U\ illustrates bounded charge-balanced minimum-power controls for spiking the same 
neuron system at various spiking times which are grater and smaller than its natural spiking 
period T = 2n. We present controls driving the neuron from 9 = to 9 = 2n at various times, 
T = 5.2, 5.3, 6.0, 7.0, 7.8, 8.2. There exist three structurally different controls which have four 
switches, two switches, and zero switches, depending on the desired spiking time. 
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Fig. 7. Optimal charge-balanced controls of minimum power given the control bound M = 0.4 for spiking a SNIPER neuron 
with w = 1 and z d = 1 at T = 5.2, 5.3, 6.0, 7.0, 7.8, 8.2. 



C. Theta Neuron Phase Model 

The theta neuron phase model is defined by f(9) = 1 + cos^ + (1 — cos 6)I b and g(Q) = 
(1 — cos 9), where I b is known as the neuron baseline current [17J. If I b > 0, then the neuron 
spikes with the period T = it / ' \[T h in the absence of any external current lit). When I b < 0, 
the neuron does not spike autonomously but it can be fired by the use of an input I(t). Since 
for I b > this neuron model can be transformed to the SNIPER phase model by a coordinate 
transformation ifTTl . we focus here on the case of < 0. Similarly, the unbounded and bounded 
charge-balanced minimum-power controls can be directly calculated by employing d8]), ©, and 
(TTOl) . or (fl4l) . (fT3T) . and (fT6l) in Section HH respectively. Optimal controls for spiking a theta 
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neuron with I b = —0.25 and M = 1 are shown in Fig. [8j 




Fig. 8. Optimal charge-balanced controls with bound M = 1.0 and for Theta neuron model (with /(, = —0.25) to elicit spikes 
at T — 4.7, 6.0,7.5, 10.0. 



The above phase models, though commonly used, are ideal mathematical models of neuron 
oscillators. We now apply our optimal control strategies to models with experimentally observed 
PRC's, such as Hodgkin-Huxley and Morris-Lecar phase models, to demonstrate their applica- 
bility and generality. 



D. Morris-Lecar Phase Model 

The Morris-Lecar model was originally proposed to capture the oscillating voltage behavior 
of giant barnacle muscle fibers (see Appendix |B]) 0. Over the past years this model has been 
extensively studied and used as a standard model for representing many different real neurons 
that are experimentally observable. For example, it has been found that Morris-Lecar PRC is 
extremely similar to the experimentally observed PRC's of Aplysia motoneuron E3l . The phase 
model of the Morris-Lecar neuron is given by 

d = u + Z(0)I(t), (26) 

where cu is the natural oscillation frequency and Z{9) represents the PRC which can be calculated 
numerically from the ODE system in Appendix |B] by the software package XPP [|24ll . For the 
set of parameter values given in Appendix [Bj the natural frequency uml = 0.283 rad/ms and 



the PRC is depicted in Fig. |9(a)| The charge-balanced minimum-power controls that elicit spikes 
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for this phase model at various times are shown in Fig. |9(b)| We consider six different cases for 
which the optimal controls have zero, two, and four switchings for spiking times that are longer 
and shorter than the natural spiking time, T = 2tt/ujml = 22.202 ms. 
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Fig. 9. |(a)| The Morris Lecar PRC for the parameters given in 
power for spiking a Morris-Lecar neuron at T = 20.5, 20.7, 21 
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(b) 

Appendix [B] |(b)| Optimal charge-balanced controls of minimum 
.0, 23.5, 24.1, 24.3 ms given the control bound M = 0.01 fiA. 



E. Hodgkin-Huxley Phase Model 

The Hodgkin-Huxley neuron model is a four dimensional system that describes the propagation 
and initiation of the action potential in squid axon (see Appendix |C]) 0. The phase model for 
this neuron oscillator is also of the form as in (|26l) . For the set of parameter values given in 
Appendix O the system has a natural frequency uhh = 0.4292 rad/ms and its PRC is displayed 



in Fig. 10(a) The charge-balanced minimum-power controls that elicit spikes at different time 



instances are shown in Fig. 10(b) 



Finally, we verified these optimal controls derived with the maximum principle by using the 
Legendre pseudospectral method. This computational method is a direct and powerful method for 
solving continuous -time optimal control problems. The basic principle is described in Appendix 
iDl and those readers interested in this method can refer to the recent comprehensive work in 
this area ll25l . Il26l . Il27ll . 11281 . The optimal controls generated by this pseudospectral method 
are presented in Fig{TH which show excellent agreement with the theoretically calculated ones 



given in Fig j 10(b) 



Phase models characterize the reduced dynamic behavior of the underlying oscillating systems, 
where the phase, but not the full state, can be observed. There is a fundamental need to explore 
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Fig. 10. 
minimum 
M = 1.0 



|(a) | The Hodgkin Huxley PRC for the parameters given in Appendix [C] |(b)| Optimal charge-balanced controls of 
power for spiking a Hodgkin-Huxley neuron at T = 13.2, 13.5, 14.0, 16.0, 16.5, 16.9 ma given the control bound 
mA. 
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Fig. 11. |(a)| Charge-balanced minimum-power controls generated by the Legendre pseudospectral method, which show excellent 
agreement with the theoretically calculated optimal controls as shown in Fig. |10(b)| |(b)| The corresponding optimal phase 
trajectories for the Hodgkin-Huxley neuron. 



the limits of the phase-reduced model as an approximation to the original oscillating system, 
because this important validation is largely lacking in the literature. The optimal controls for 
phase models presented so far in this work change the spiking times of an oscillator during the 
course of one oscillatory cycle, so that a desired spike train can be constructed by repeating the 
control input. We now apply the optimal controls derived according to the scalar Hodgkin-Huxley 
phase model to its full state-space model, which is a system of four differential equations as 
shown in Appendix The spike train obtained by repeated application of the optimal control 
producing an inter-spike time T = 16 ms, subject to the control amplitude bound M = 1 mA, 
and the uncontrolled train spiking at the natural period, T = 2h/uj H h = 14.64 ms, are illustrated 
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Fig. 12. Uncontrolled and controlled spiking trains of Hodgkin-Huxley Model 



in Fig. [[21 It is seen that the optimal control delays the spiking time from 14.64 ms to 16.02 ms 
in the state-space model. 

IV. Discussion 

In this paper, we considered the optimal control of phase models of neuron oscillators. We 
derived charge-balanced minimum-power current stimuli that elicit spikes of neurons at desired 
time instances for the cases of unbounded and bounded control amplitude. In particular, we 
showed that for the bounded case the optimal control has switching characteristics synthesized 
by the unbounded optimal control and the control bound. We implemented the resulting analytical 
optimal controls to various commonly used phase models, including mathematically ideal and 
experimentally observed models, to demonstrate their applicability. We then applied the optimal 
controls derived according to the phase-reduced model of Morris-Lecar and Hodgkin-Huxley to 
the corresponding full state-space system to validate the approximation of the phase model under 
weak forcing. The theory presented in this work can be applied not only to neuron oscillators but 
also to any oscillating systems that can be represented using similar model reduction techniques 
such as biological, chemical, electrical, and mechanical oscillators. 

The theoretical results presented in this paper characterize the fundamental limit of how the 
dynamics of neurons can be perturbed by the use of external inputs. Alternatively, they provide 
an insight into how the neuron dynamics determine the synaptic input necessary for eliciting 
spikes, which facilitates the development of optimal stimuli for neurological treatments such as 
deep brain stimulation for Parkinson's disease. The extension of this work to the optimal control 
of networks of neuron oscillators is of fundamental and practical importance. Our recent work 
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has shown that a ensemble of uncoupled neurons is controllable, and that the minimum-power 
controls that spike a network of heterogeneous neurons can be found by using a multidimensional 
pseudo spectral method |fT9ll . We plan to extend this recent work to investigate controllability of 
coupled neurons and related optimal control problems. Systems described by the Kuramoto 
model will be considered. 



Appendix A 

Spiking Sinusoidal Neurons with Bounded Control 



Simple first and second order optimality conditions applied to (|T8l) find that the maximum 



value of I* occurs at 6 = tt/2 for c < and at 6 = 3n/2 for c > (see Fig ]13(a)| for c < 0). 
According to (fl9l) . c = corresponds to T = 2tv/uj and c < (c > 0) corresponds to T < 2n/uj 
(T > 2tx I uS). Therefore, the constant c for the shortest spiking time with the control I* satisfying 
\I*(t)\ < M can be calculated by substituting I* = M and 9 = n/2 to £[8]), and then from (03 
we obtain the shortest spiking period by J*, T m \ n , as in (1231) . 



nil n 3n/2 Itc 6 nil B n 9 371/2 9 In 

Phase (6) 1 2 Phase (e) 3 4 

(a) (b) 

Fig. 13. |(a)| An illustration of the optimal control J* with its maximum value occurring at 6 = 7r/2 for c > 0, which gives the 
shortest possible spiking time subject to the control bound M. |(b)| An illustration of the case when I* > M with intersections 
at 0i, 2 , 03, and 4 . 



Since I* takes the maximum value at 6 = 3n/2 for c > 0, which corresponds to T > 2ti/uj, 
we have |/*| < (u — v/a> 2 — czfy/zd, which leads to |/*| < uo / Zd < M for T > 2ti/lu, provided 
that M > uj/Zd- This implies that I* is the minimum-power control for any desired spiking time 
T > 2n/cu when M > to/zd. Since the smallest spiking time by the control I* with the bound M 
is given by T m \ n as in (|23l) . I* as described in (TT8l) and (|T9l ) will be the optimal control for any 
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spiking time T > T^ in if the bound satisfies M > u/z d . A shorter spiking time T G [T^ in , T^ in ) 
is feasible but can not be achieved by I* alone. Suppose that T G [T^f in ,T^ in ), then there exist 
two angles 9\ = sin _1 [— 2Mu/ (z d M 2 + cz d )] and 9 2 = n — 9i where I* meets the bound M, 
illustrated in Fig. 13(b). When 9 G (9 t ,9 2 ), I* > M and we take 1(9) = M for 9 G [9x,9 2 ], 
Then, from (fl"3T) . the Lagrange multiplier is A = (H — M 2 )/(u + z d M sin 9). This multiplier 
satisfies the adjoint equation ©, therefore 1(9) = M is optimal for 9 G [#i,#2]- Similarly, by 
symmetry, I* < — M when 9 G [#3, #4], where 3 = n + 6>i and 9 4 = 2n — 9\, if the desired 
spiking time is T G [T^f in: T^ in ). It can be easily shown by the same fashion that 1(9) = —M 
is optimal in the interval 9 G [6*3, 6*4]. Therefore, the minimum-power optimal control that spikes 
the neuron at T G [T^{ n , T^* n ) can be characterized by four switchings between I* and M as 
shown in (|24l) . 

Appendix B 
Morris-Lecar Model 

The dynamics of the Morris-Lecar neuron are described by two coupled dynamical equations 

V = i [(I b + I)+ gc a moo(Vca ~V)+ g k w(V k ~V) + g L (V L - V)} 

W = ^(Woo - w)/t w (V) 

moo = 0.5[1 + tanh((V - I4)/V 2 )] 
Woo = 0.5[1 + tanh((y - V 3 )/V 4 )) 
T u = l/co8h[(V-V 3 )/(2Vi)]. 
In Section IIII-Dl we consider the following parameter values 

(j) = 0.5, I b = 0.09 LiA/cm 2 , V x = -0.01 mV 

V 2 = 0.15 mV, V 3 = 0.1 mV, V 4 = 0.145 mV, 

gca = 1 mS/cm 2 , Vk = —0.7 mV, Vl = —0.5 mV, 
gk = 2 mS/cm 2 , gL = 0.5 mS/cm 2 , C = 1 fiF/cm 2 
Vca = 1 mV^ 



September 20, 2011 



DRAFT 



20 



Appendix C 

The dynamics of the Hodgkin-Huxley neuron are described by a set of differential equations 

CV — I — -g Na h(V - V Na )m 3 - g k (V - V k )n 4 - g L (V - V L ) 
m = a m (V)(l — m) — b m (V)m 
h = a h (V)(l - h) - b h (V)h 
h = a n (V)(l - n) - b n (V)n 
a m (V) = 0.1(V + 40)/[l - exp(-(V + 40)/10)] 
6 m (y) = 4e X p[-(y + 65)/18] 
a h (V) = 0.07exp[-(V + 65)/20] 
6 h (V) = l/(l + exp[-(y + 35)/10)] 
a n (y) = 0.01(V + 55)/[l - exp(-(V r + 55)/10)] 
b n (V) = 0.125 exp[-(y + 65)/80]. 

In Section IIII-EL we consider the following parameter values 

V Na = 50 mV, V k = -77 mV, v L = -54.4 mV, 

gNa = 120 mS/cm 2 , g k = 36 mS/cm 2 , gi = 0.3 mS/cm 2 , 
C = 1 fiF/cm 2 , I = 10 fiA/cm 2 . 

Appendix D 

Legendre pseudospectral method for optimal control of phase-reduced 

oscillators 

The pseudospectral method is a spectral collocation method that was originally developed 
to solve partial differential equations, and has recently been adapted to solve optimal control 
problems H2^,E9,lE2,E8,lE3,|23.Inthis approach the differential equations that relate 
the states and the controls are discretized at specific collocation nodes, which results in a discrete 
optimization problem. All continuous-time functions are rescaled to the time domain of [-1,1] 
and expanded by an orthogonal polynomial basis based on a set of selected quadrature nodes 
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||29l . Here, we use the Legendre-Gauss-Lobatto(LGL) nodes, and can then write the iVth order 
interpolating approximations of the state and control functions 



N 



6(t)aI N 9(t) = ^0 fc 4(t), 

k=0 

N 

I(t)*I N I(t) =^44(t), 



k=0 



where 



N 



.•(*)= n 



t-U 



k = 0,l,..-,N, 



i=0,i^k 

are the Lagrange polynomials with 4(£j) — hi, the Kronecker delta function. The derivative of 
lN0{t) at the LGL node tj, j = 0, 1, 2, . . . , N is then given by 

. N N 

a 



—Is0{tj) — ^^6kik(tj) — DjkOk, 



k=0 



k=0 



where Dj k are the jfc t/l elements of the constant (N + 1) X (N + 1) differentiation matrix defined 
by 

L N(tj) 1 • / I 

L N (t k ) tj-t k J I 
— JV^JV+1) J = k = Q 

j = k = N 
otherwise. 



D 



N(N+1) 
4 



The integral cost functional of the optimal control problem as in © can be accurately ap- 
proximated by the Gauss-Lobatto integration rule. Thus, the pseudospectral discretization of the 
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optimal control problem © gives rise to a nonlinear program of the form 



min - I? w i, 
i ...i n I . =0 

N T 

s.t. Y,D ik e k = -[f(e i ) + i i g(e i )] 

k=0 

N rp 

D ikPk = -^h, 



k=0 








e = 


0. 


On 


= 2tt 


Po = 


o. 


Pn 


= 0, 



\h\ < M, 

where Wi are the LGL weights given by Wi = N ^+i) {L N \t )) 2 - Solvers for this type of mini- 
mization problems are readily available and straightforward to implement. We approximate the 
problem using 151 nodes (N = 150) and implement it in the AMPL language [|30l . We use a 
third party nonlinear programming solver KNITRO [31 J to solve this optimization. This Legendre 
pseudospectral method provides a direct method to verify the analytical results presented in 
Section HI 
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